In [6]:
from __future__ import division
import time
import numpy as np
from pylab import ion
import matplotlib as mpl
import matplotlib.pylab as p
from matplotlib import pyplot as plt
from matplotlib import animation
#from numba import jit, autojit
ion()
Define the successive overrelaxation algorithm's "step" (each gridpoint gets updated once)
In [7]:
#@autojit
def SOR_step(vx, vy, P, Nx, Ny, boundaries=[], omega=1):
#for i in range(1,Nx-1):
for i in range(Nx-2,0,-1):
#for j in range(1,Ny-1):
for j in range(Ny-2,0,-1):
if (i,j) not in boundaries:
new_vx = 0.25*(vx[i+1,j] + vx[i-1,j] + vx[i,j+1] + vx[i,j-1] \
-h/2.*vx[i,j]*(vx[i+1,j]-vx[i-1,j]) \
-h/2.*vy[i,j]*(vx[i,j+1]-vx[i,j-1]) \
-h/2.*(P[i+1,j]-P[i-1,j])
)
residual = new_vx - vx[i,j]
vx[i,j] = vx[i,j]+omega*residual
Parameters
In [8]:
Nx = 200
Ny = 200
h = 1
L = 3
nu = 1
rho = 10e3
V0 = 1
domain_width = 10 # m
domain_height = 2 # m
omega = 1.2
Set initial velocities and pressures in entire region
In [9]:
vx = np.zeros((Nx, Ny))
vy = np.zeros((Nx, Ny))
P = np.zeros((Nx,Ny))
In [10]:
pixel_width = domain_width/Nx
pixel_height = domain_height/Ny
plate_y0 = int(np.round(Ny/2+h/2/pixel_height))
plate_y1 = int(np.round(Ny/2-h/2/pixel_height))
plate_x0 = int(np.round(Nx/2-L/2/pixel_width))
plate_x1 = int(np.round(Nx/2+L/2/pixel_width))
Specify boundary conditions around edges
In [11]:
vx[:,0] = 1 #bottom
vx[:,-1] = 1 #top
vx[0,:] = V0 #left
vx[-1,:] = 1 #rigfht
vy[0,:] = 0
vy[-1,:] = 0
vy[:,0] = 0
vy[:,-1] = 0
P[0,:] = 0
P[-1,:] = 0
P[:,0] = 0
P[:,-1] = 0
Specify boundary conditions on plates
In [12]:
vx[plate_x0:plate_x1+1,plate_y0] = 0
vx[plate_x0:plate_x1+1,plate_y1] = 0
boundaries = []
for i in range(plate_x0, plate_x1+1):
for j in [plate_y0, plate_y1]:
boundaries.append( (i,j) )
In [13]:
#%pylab inline
Plot initial conditions
In [14]:
f1=plt.figure(1, figsize=(15,3), dpi=60)
f1.clf()
ax1=f1.add_subplot(111)
ax1.imshow(vx.T, interpolation='nearest', cmap=plt.cm.jet, origin='lower')
title(r"Intitial conditions, $v_x$")
Out[14]:
In [15]:
iteration = 0
error = 100
old_vx = vx.copy()
while np.max(error) > 0.01:
SOR_step(vx, vy, P, Nx, Ny, boundaries=boundaries, omega=omega)
error = np.abs(vx-old_vx)
old_vx = vx.copy()
# if iteration % 10 == 0:
# ax1.imshow(vx.T, interpolation='nearest', cmap=plt.cm.jet,
# vmin=0, vmax=2, origin='lower')
# plt.draw()
# plt.show()
# time.sleep(0.2)
iteration += 1
print iteration
Plot results
In [16]:
f2=plt.figure(2, figsize=(8,6.5), dpi=60)
f2.clf()
ax2=f2.add_subplot(111)
im2 = ax2.imshow(vx.T, interpolation='nearest', cmap=plt.cm.jet, origin='lower')
f2.colorbar(im2, ax=ax2)
title(r"Converged solution, $v_x(x,y)$")
xlabel(r"$x$"); ylabel(r"$y$")
tight_layout()
In [16]: